function [triangleIndex] = searchTriangle(p,sites,triangles, testIndex,gridStructure)

triangleIndex = 0;

% try testIndex first to save time
a = [sites(triangles(testIndex,1),[1 2]) 0]; % x, y axis of site, where z = 0;
b = [sites(triangles(testIndex,2),[1 2]) 0];
c = [sites(triangles(testIndex,3),[1 2]) 0];

if(sameSide(p,a,b,c) && sameSide(p,b,a,c) && sameSide(p,c,a,b))
    triangleIndex = testIndex;
    return;
end 

% if not testIndex search all triangles
Nx = 30; Ny = 30;
minx = min(sites(:,1)); maxx = max(sites(:,1));
miny = min(sites(:,2)); maxy = max(sites(:,2));
dx = (maxx-minx+1)/Nx;
dy = (maxy-miny+1)/Ny;

indx = floor((p(1) - minx)/dx) + 1;
indy = floor((p(2) - miny)/dy) + 1;

%boundary case
if(indx>Nx) indx = Nx;
elseif(indx<1) indx = 1;
end
if(indy>Ny) indy = Ny;
elseif(indy<1) indy = 1;
end

for i=1:size(gridStructure(indx,indy).triangles,1)
    testIndex = gridStructure(indx,indy).triangles(i);
    a = [sites(triangles(testIndex,1),[1 2]) 0]; % x, y axis of site, where z = 0;
    b = [sites(triangles(testIndex,2),[1 2]) 0];
    c = [sites(triangles(testIndex,3),[1 2]) 0];

    if(sameSide(p,a,b,c) && sameSide(p,b,a,c) && sameSide(p,c,a,b))
        triangleIndex = testIndex;
        return;
    end 
end

if(triangleIndex==0)
    disp(['Looking for : ' num2str(p)]);
    for i=1:size(triangles,1)
        a = [sites(triangles(i,1),[1 2]) 0]; % x, y axis of site, where z = 0;
        b = [sites(triangles(i,2),[1 2]) 0];
        c = [sites(triangles(i,3),[1 2]) 0];


        if(sameSide(p,a,b,c) && sameSide(p,b,a,c) && sameSide(p,c,a,b))
            triangleIndex = i;
            disp(['triangle found: ' num2str(triangleIndex)]);
            break;
        end 
    end
    
end

end

%%
function result = sameSide(p1,p2,a,b)
    cp1 = cross(b-a, p1-a);
    cp2 = cross(b-a, p2-a);
    
    result = dot(cp1,cp2) >= 0;
    
end

